module CImod
	use compute
	implicit none

	contains
	
	function getLR2s(y1,y2,Zsum,fa1,fa2,tst) result(val)
		real	:: y1(1:k),y2(1:k),Zsum,fa1,fa2,val
		type(tsttype)	:: tst

		type(ptype)	:: p
		real		:: h
		integer		:: i
					
		call setpair([Zsum,Y1],[0.0,Y2],p)
		h=0
		do i=1,tst%n0active
			h=h+tst%lam(i)*getdens0p_stage2(tst%thlist0(:,i),p)
		enddo
		if(h/(fa1*fa2)<1E-20) then
			val=-100
		else
			val=log(h/(fa1*fa2))
		endif
	end function
	
	function getLR1s(y1,y2,Zsum,fa1,tst) result(val)
		real	:: y1(1:k),y2(1:k),Zsum,fa1,val
		type(tsttype)	:: tst
		real		:: z,h,ma,msx(3),var
		integer		:: j
					
		Z=Zsum-sum(y2)
		var=1+sum(y2**2)
		h=0
		do j=1,tst%n0active_stage1
			msx=tst%thlist0_stage1(:,j)
			ma=getmeanadj(msx,y1(k))
			h=h+tst%lam_stage1(j)*getdens_base(msx,y1)*exp(-0.5*(z+ma)**2/var)/sqrt(var)
		enddo
		if(h/fa1<1E-20) then
			val=-100
		else
			val=log(h/fa1)
		endif
		val=val-5.0*switchind(y2)
	end function
	
	function gett2(y1,y2,Zsum,tst) result(val)
		real	:: y1(1:k),y2(1:k),Zsum,val
		type(tsttype)	:: tst
		real	:: var,cv,cvw

		var=1.0+sum(y1**2)+sum(y2**2)
		cvw=1.0/var
		cv=cvw*tst%cv_zt(1)+(1.0-cvw)*tst%cv_zt(2)
		val=cv**2-(Zsum+sum(y1)-sum(y2))**2/var
	end function		
		
	function getLRs(y1,y2,Zsum,tst) result(val)
		real	:: y1(k),y2(k),Zsum,val(4)
		type(tsttype)	:: tst
		real	:: fa1,fa2
		fa1=getf1(y1); fa2=getf1(y2)

		val(1)=getLR2s(y1,y2,Zsum,fa1,fa2,tst)
		val(2)=getLR1s(y1,y2,Zsum,fa1,tst)
		val(3)=getLR1s(y2,y1,-Zsum,fa2,tst)
		val(4)=gett2(y1,y2,Zsum,tst)
	end function
	
	function getcross(x,y) result(val)
		real	:: x(3),y(3),val
		real	:: a,b,c,d,xp
		val=1000
		d=((x(1) - x(2))*(x(1) - x(3))*(x(2) - x(3)))
		a=(x(3)*(-y(1) + y(2)) + x(2)*(y(1) - y(3)) + x(1)*(-y(2) + y(3)))/d
		b=(x(3)**2*(y(1) - y(2)) + x(1)**2*(y(2) - y(3)) + x(2)**2*(-y(1) + y(3)))/d
		c=(x(3)*(x(2)*(x(2) - x(3))*y(1) + x(1)*(-x(1) + x(3))*y(2)) + x(1)*(x(1) - x(2))*x(2)*y(3))/d
		d=b**2-4*a*c
		if(d<0) return
		val=(-b+sqrt(d))/(2*a)
		if(val>=x(1)-1E-6 .and. val<=x(3)+1E-6) return
		val=(-b-sqrt(d))/(2*a)
		if(val>=x(1)-1E-6 .and. val<=x(3)+1E-6) return
		val=1000
	end function
	
	function getCIlb(y1,y2,Zsum,tst) result(val)
		real	:: y1(k),y2(k),Zsum,val
		type(tsttype)	:: tst
		real	:: mu,s,cmus(3),cLRs(4,3),maxLR
		real	:: mustep
		integer	:: i,j
		
		mu=(sum(Y1)-sum(Y2)+Zsum)/n_sample
		s=sqrt(1+sum((Y1-mu)**2)+sum((Y2+mu)**2))/real(n_sample)
		mustep=max(0.05*s,0.5/(n_sample-2*k))
		mu=mu-6*s
		do
			mu=mu-2*s
			cmus(2)=mu
			cLRs(:,2)=getLRs(y1-mu,y2+mu,Zsum-(n_sample-2*k)*mu,tst)
			if(maxval(cLRs(:,2))<-5.0) exit 
		enddo
		mu=mu+mustep
		cmus(3)=mu
		cLRs(:,3)=getLRs(y1-mu,y2+mu,Zsum-(n_sample-2*k)*mu,tst)
		do j=1,500
			maxLR=maxval(cLRs(:,3))
			mu=mu+merge(.5,merge(1.0,2.0,maxLR>-4.0),maxLR>-2.0)*mustep
			cmus=[cmus(2:3),mu]
			cLRs=reshape([cLRs(:,2:3),getLRs(y1-mu,y2+mu,Zsum-(n_sample-2*k)*mu,tst)],[4,3])
			val=1000
			do i=1,4
				val=min(val,getcross(cmus,cLRs(i,:)))
			enddo

			if(val<1000) return
		enddo
		print *,"couldn't find CI lower bound"
		call mdisp([y1,y2,Zsum])
		val=(sum(Y1)-sum(Y2)+Zsum)/n_sample
		return
	end function
	
end module